home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
c
/
fftsing.zip
/
SING.C
< prev
next >
Wrap
C/C++ Source or Header
|
1992-04-15
|
18KB
|
844 lines
/**************************************************************************
Javier Soley, Ph. D, FJSOLEY @UCRVM2.BITNET
Escuela de Física y Centro de Investigaciones Geofísicas
Universidad de Costa Rica
***************************************************************************/
/* Computes the DISCRETE FOURIER TRANSFORM of very long data series.
* Restriction: the data has to fit in conventional memory.
*
* Compile in compact or large models ===> Pointers must be FAR
* Two huge pointers are used to access the real and imaginary
* parts of the transform without 64k wrap around.
*
* This functions are translations from the fortran program in
*
* R. C. Singleton, An algorithm for computing the mixed radix fast
* Fourier transform
*
* IEEE Trans. Audio Electroacoust., vol. AU-17, pp. 93-10, June 1969.
* Some features are:
*
* 1-) Accepts an order of transform that can be factored not only
* in prime factors such 2 and 4, but also including odd factors
* as 3, 5, 7, 11, etc.
* 2-) Generates sines and cosines recursively and includes
* corrections for truncation errors.
* 3-) The original subroutine accepts multivariate data. This
* translation does not implement that option (because I
* do not needed right now).
*
* Singleton wrote his subroutine in Fortran and in such a way that it
* could be ported allmost directly to assembly language. I transcribed
* it to C with little effort to make it structured. So I apologize to
* all those C purists out there!!!!!!!!
*
*/
/* Version 2.0 March/30/92 */
/* Includes */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* Defines */
#define TWO_PI ((double)2.0 * M_PI)
#define MAXF 23
#define MAXP 209
#define MY_SIZE double /* or float or long double */
#define MY_TYPE huge /* or far if data fits in two 64 K segments */
/* Globals */
long nn, m, flag, jf, jc, kspan, ks, kt, nt, kk, i;
MY_SIZE c72, s72, s120, cd, sd, rad, radf, at[23], bt[23];
long nfac[23];
int inc;
long np[MAXP];
/* The functions */
void radix_2( a, b)
MY_SIZE MY_TYPE *a;
MY_SIZE MY_TYPE *b;
{ long k1, k2;
MY_SIZE ak, bk, c1, s1;
kspan >>= 1;
k1 = kspan +2;
do {
do {
k2 = kk + kspan;
ak = a[k2-1];
bk = b[k2-1];
a[k2-1] = a[kk-1] -ak;
b[k2-1] = b[kk-1] -bk;
a[kk-1] += ak;
b[kk-1] += bk;
kk = k2 + kspan;
} while ( kk <= nn);
kk = kk - nn;
} while ( kk <= jc);
if ( kk > kspan) flag = 1;
else
{
do {
c1 = 1.0 - cd;
s1 = sd;
do {
do {
do {
k2 = kk + kspan;
ak = a[kk-1]- a[k2-1];
bk = b[kk-1]- b[k2-1];
a[kk-1] += a[k2-1];
b[kk-1] += b[k2-1];
a[k2-1] = c1*ak - s1*bk;
b[k2-1] = s1*ak + c1*bk;
kk = k2 + kspan;
} while ( kk < nt );
k2 = kk - nt;
c1 = -c1;
kk = k1 - k2;
} while ( kk > k2 );
ak = c1- (cd*c1+sd*s1);
s1 = (sd*c1-cd*s1) +s1;
/***** Compensate for truncation errors *****/
c1 = 0.5/(ak*ak+s1*s1)+0.5;
s1 *= c1;
c1 *= ak;
kk += jc;
} while ( kk < k2);
k1 = k1 + inc + inc;
kk = (k1- kspan) /2 + jc;
} while ( kk <= jc + jc );
}
}
void radix_4( isn, a, b)
MY_SIZE MY_TYPE *a;
MY_SIZE MY_TYPE *b;
int isn;
{
long k1, k2, k3;
MY_SIZE akp, akm, ajm, ajp, bkm, bkp, bjm, bjp;
MY_SIZE c1, s1, c2, s2, c3, s3;
kspan /= 4;
cuatro_1:
c1 = 1.0;
s1 = 0;
do {
do {
do {
k1 = kk + kspan;
k2 = k1 + kspan;
k3 = k2 + kspan;
akp = a[kk-1] + a[k2-1];
akm = a[kk-1] - a[k2-1];
ajp = a[ k1-1] + a[k3-1];
ajm = a[ k1-1] - a[k3-1];
a[kk-1] = akp + ajp;
ajp = akp - ajp;
bkp = b[kk-1] + b[k2-1];
bkm = b[kk-1] - b[k2-1];
bjp = b[k1-1] + b[k3-1];
bjm = b[k1-1] - b[k3-1];
b[kk-1] = bkp + bjp;
bjp = bkp - bjp;
if ( isn < 0) goto cuatro_5;
akp = akm - bjm;
akm = akm + bjm;
bkp = bkm + ajm;
bkm = bkm - ajm;
if (s1 == 0.0) goto cuatro_6;
cuatro_3: a[ k1-1] = akp*c1 - bkp*s1;
b[ k1-1] = akp*s1 + bkp*c1;
a[ k2-1] = ajp*c2 - bjp*s2;
b[ k2-1] = ajp*s2 + bjp*c2;
a[ k3-1] = akm*c3 - bkm*s3;
b[ k3-1] = akm*s3 + bkm*c3;
kk = k3 + kspan;
} while ( kk <= nt);
cuatro_4:
c2 = c1 - (cd*c1 + sd*s1);
s1 = (sd*c1 - cd*s1) + s1;
/***** Compensate for truncation errors *****/
c1 = 0.5 / (c2*c2 + s1*s1) +0.5;
s1 = c1 * s1;
c1 = c1 * c2;
c2 = c1*c1 - s1*s1;
s2 = 2.0 * c1 *s1;
c3 = c2*c1 - s2*s1;
s3 = c2*s1 + s2*c1;
kk = kk -nt + jc;
} while ( kk <= kspan);
kk = kk - kspan + inc;
if ( kk <= jc) goto cuatro_1;
if ( kspan == jc) flag =1;
goto out;
cuatro_5:
akp = akm + bjm;
akm = akm - bjm;
bkp = bkm - ajm;
bkm = bkm + ajm;
if (s1 != 0.0) goto cuatro_3;
cuatro_6:
a[k1-1] = akp;
b[k1-1] = bkp;
b[k2-1] = bjp;
a[k2-1] = ajp;
a[k3-1] = akm;
b[k3-1] = bkm;
kk = k3 + kspan;
} while ( kk <= nt);
goto cuatro_4;
out:
s1 = s1 + 0.0;
}
/* Find prime factors of n */
void fac_des( n )
long n;
{
long k, j, jj, l;
k = n;
m = 0;
while ( k-(k / 16)*16 == 0 ) {
m++;
nfac[m-1] = 4;
k /= 16;
}
j = 3;
jj = 9;
do {
while (k % jj == 0) {
m++;
nfac[m-1] = j;
k /= jj;
}
j += 2;
jj = j * j;
} while ( jj <= k);
if (k <= 4) {
kt = m;
nfac[m] = k;
if (k != 1) m++;
}
else {
if (k-(k / 4)*4 == 0) {
m++;
nfac[m-1] = 2;
k /= 4;
}
kt = m;
j = 2;
do {
if (k % j == 0 ) {
m++;
nfac[m-1] = j;
k /= j;
}
j = ((j+1)/ 2)*2 + 1;
} while ( j <= k);
}
if (kt != 0) {
j = kt;
do {
m++;
nfac[m-1] = nfac[j-1];
j--;
} while ( j != 0);
}
}
/* Permute the results to normal order */
void permute( ntot, n, a, b)
long ntot, n;
MY_SIZE MY_TYPE *a;
MY_SIZE MY_TYPE *b;
{
long k, j, k1, k2, k3, kspnn, maxf;
MY_SIZE ak, bk;
long ii, jj;
maxf = MAXF;
np[0] = ks;
if (kt != 0) {
k = kt +kt +1;
if (m < k) k--;
j = 1;
np[k] = jc;
do {
np[j] = np[j-1] / nfac[j-1];
np[k-1] = np[k] * nfac[j-1];
j++;
k--;
} while (j < k);
k3 = np[k];
kspan = np[1];
kk = jc+1;
k2 = kspan + 1;
j = 1;
/***** Permutation of one dimensional transform *****/
if (n == ntot) {
do {
do {
ak = a[kk-1];
a[kk-1] = a[k2-1];
a[k2-1] = ak;
bk = b[kk-1];
b[kk-1] = b[k2-1];
b[k2-1] = bk;
kk += inc;
k2 += kspan;
} while ( k2 < ks);
ocho_30: do {
k2 -= np[j-1];
j++;
k2 += np[j];
} while (k2 > np[j-1]);
j = 1;
ocho_40: j = j + 0;
} while (kk < k2);
kk += inc;
k2 += kspan;
if (k2 < ks) goto ocho_40;
if (kk < ks) goto ocho_30;
jc = k3;
}
else { /* Permutation for multiple transform */
ocho_50:
do {
do {
k = kk + jc;
do {
ak = a[kk-1];
a[kk-1] = a[k2-1];
a[k2-1] = ak;
bk = b[kk-1];
b[kk-1] = b[k2-1];
b[k2-1] = bk;
kk += inc;
k2 += inc;
} while ( kk < k);
kk = kk + ks - jc;
k2 = k2 + ks - jc;
} while (kk < nt);
k2 = k2 - nt + kspan;
kk = kk - nt + jc;
} while (k2 < ks);
do {
do {
k2 -= np[j-1];
j++;
k2 += np[j];
} while (k2 > np[j-1]);
j =1;
do {
if ( kk < k2 ) goto ocho_50;
kk +=jc;
k2 += kspan;
} while (k2 < ks);
} while (kk < ks);
jc = k3;
}
}
if ( (2*kt +1) < m) {
kspnn = np[kt];
/* Permutation of square-free factors of n */
j = m - kt;
nfac[j] = 1;
do {
nfac[j-1] *= nfac[j];
j--;
} while (j != kt);
kt++;
nn = nfac[kt-1] -1;
if (nn > MAXP) {
printf("product of square free factors exceeds allowed limit\n");
exit(2);
}
jj =0;
j=0;
goto nueve_06;
nueve_02:
jj -= k2;
k2 = kk;
k++;